function y = LWPR(q, X, Y, M, h)

% Locally weighted polynomial regression
% q, query point
% X, input points
% Y, output points
% M, order of polynomials used for regression
% h, width of gaussian kernel, larger means giving (near) equal weight to all points, i.e. globally weighted regression

% setup variables
WX = [ones(size(X, 1), 1) X];
WY = Y;
q = [1 q];

% compute weights, scale inputs and outputs
ws = sqrt(exp(-dist2(q, WX)/h));
for i = 1:length(ws)
    WX(i, :) = WX(i, :) * ws(i);
    WY(i, :) = WY(i, :) * ws(i);
end

% compute polynomial
polynomial = computePolynomial(M, size(WX, 2));
TWX = computeTerms(WX, polynomial);
q = computeTerms(q, polynomial);

% compute coefficients
B = inv(TWX' * TWX) * TWX' * WY;

% perform regression
y = q * B;
